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We describe a simulation method for a quantum spin model of a generic, general purpose 
quantum computer. The use of this quantum computer simulator is illustrated through several 
implementations of Grover's database search algorithm. Some preliminary results on the stability 
of quantum algorithms are presented. 
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§1. Introduction 

The idea tliat a Quantum Computer (QC) might be 
more powerful than an ordinary computer is based on the 
notion that a quantum system can be in any superposi- 
tion of states and that interference of these states allows 
exponentially many computations to be done in paral- 
lehllP A QC may solve certain computationally hard prob- 
lems such as factoring integers and snfirdiiiie databases 
faster than a conventional computer .Hu qIjO Q' This in- 
trinsic parallelism might be used to solve other difficult 
problems as well, such as for example the calculation 
of thaghjjsical properties of quantum many-body sys- 
tems.EloEj^tj-' In fact, part of Feynman's original moti- 
vation to consider QC's was that they might be used as 
a vehicle to performj^act simulations of quantum me- 
chanical phenomena.^ 

Theoretical work on quantum computation usually as- 
sumes the existence of units that perform highly idealized 
unitary operations. However, in practice these opera- 
tions are difficult to realize. Disregarding decoherence, a 
hardware implementation of a QC will perform unitary 
operations that are more complicated than those con- 
sidered in most theoretical work: In a QC the internal 
quantum dynamics of each elementary constituent is a 
key ingredient of the QC itself. 

This paper describes a simulator for a generic physical 
model of a QC, strictly working according to the laws of 
quantum mechanics. We implement Grover's database 
search quantum algorithm (QA)I3'LP using ideal and more 
realisJ,kL,mii±a. ^iuch as those used in the 2-qubit NMR 



§2. Quantum Spin Dynamics 

Gcnerically, hardware QC's are modeled in terms of 
S=l/2 spins (qubits) that evolve in time according to 
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the time-dependent Schrodinger equation (TDSE) 

^^^\m)-Hit)\m), (2.1) 

in units such that h = 1 and where describes the 

state of the whole QC at time t. The time-dependent 
Hamiltonian H{t) takes the form 



H{t) 



E 



E 



Jj;k,a{t)S°'S'j^ 



j,k=l a=x,y,z 



-EE if^i^-At) 

j = l a=x,y,z 

sin(/j- „i + Sf, (2.2) 

where the first sum runs over all pairs P of spins, 5" de- 
notes the a-th component of the spin-1/2 operator repre- 
senting the j-th qubit, Jj,k,a{t) determines the strength 
of the interaction between the qubits labeled j and fc, 
hj,a,oit) and hj^a,i(t) are the static (magnetic) and peri- 
odic (RF) field acting on the j-th spin respectively. The 
frequency and phase of the periodic field are denoted by 
fj,a and Lpj,a- The number of qubits is L and the di- 
mension of the Hilbert space D ~ 2^. Hamiltonian (2.2) 
captures the physics of most candidate technologies for 
building QC's. 

A QA for the QC modeled by ( [2I] ) and (2.2) consists 
of a sequence of elementary operations (EO's). The ac- 
tion of an EO on the state 1 5*) of the quantum processor 
is determined by the values of all the J's and /I's (which 
are kept constant during the operation) and the time in- 
terval it is active. The input state |^'(i)) is transformed 
into the output state |5'(t -I- r)) where r denotes the time 
it takes to complete the EO. During this time interval the 
only time-dependence of H{t) is through the (sinusoidal) 
modulation of the fields on the spins. The time evolution 
of the QC itself is governed by the TDSE (O). 
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Formally the solution of ( p.lD can be expressed in 
terms of the unitary transformation U{t + T,t) = 



exp_|_(— ? H{u)du), where exp^ denotes the time- 
ordered exponential function. Using the semi-group 
property of U{t + r, t) we can write 



U{t + T,t) = U{t + m5,t+{m- 1)5) ■ ■ 
U{t + 2S,t + S)U{t + S,t), 



(2.3) 



where r = mS (m > 1). The standard procedure to 



construct algorithms for solving the TDSE (2.1) is to 
replace each U{t + (n + 1)6, t -t^jpAIJay a symmetrized 
product-formula approximationlla^tSEH) For the casBySd 
hand a convenient choice is (other decomposition^' £3 
work equally well but are somewhat less efficient for our 
purposes): 

U{t +{n + 1)6, t + n6)^ U{t + {n + 1)6, t + n6), (2.4) 



U(t+ {n + l)6,t + n6) = 

^-iSH,{t+{n+l/2)S)/2^-iSHy(t+(n+l/2)S)/2 
^^-iSH^{t+{7i+l/2)S)^-iSHy{t+(n+l/2)S}/2 



xe 



-iSH^ (t+(n+l/2)S) /2 



(2.5) 



where 



H^it) 

L 

-E 

j.k=i 



[hj.afi + hj,a,l Sin{fj,at + ipj^a)) Sf, (2.6) 



with a — x,y, z. In (2.6) the time dependence of the J's 
and the /I's has been omitted because these parameters 
are constant during the execution of an EO. 

Evidently \J(t ~\- t, t) is unitary by construction, im- 
plying that thej-algorithm to solve the TDSE is uncondi- 
tionally stableJlS^ It is easily shown that the-.algorithm 
is correct to second order in the time-step 6x3 Further- 
more U{t + T, t) can be used aa-a biulding block to con- 
struct higher-order algorithms.aESEj'H In practice it is 
easy to find reasonable, relatively small, values of m such 
that the results obtained no longer depend on m (and 6). 
Then, for all practical purposes, these results arc indis- 
tinguishable from the exact solution of the TDSE ( ^.1[ ) . 

It is customary to take as basis states {I'/'n)} the direct 
product of the eigenvectors of the 5| (i.e. spin-up 
and spin-down \l)j). In this basis e-*'5^'(*+("+i/2)'5)/2 
changes the input state by altering the phase of each 
of the basis vectors. As is a sum of pair interac- 
tions it is trivial to rewrite this operation as a direct 
product of 4x4 diagonal matrices (equivalent to the so- 
called interaction-controlled phase shifts) and 4x4 unit 
matrices. Hence the computation of eyip{—i6Hz{t -I- (n -|- 
l/2)(5)/2)|^') has been reduced to the multiplication of 
two vectors, element- by-element. The unitary matrix 
^-iSHy(t+(n+i/2)S)/2 written in a similar manner 



but the matrices that contain the interaction-controlled 
phase-shift have to be replaced by non-diagonal matri- 
ces. Although this does not present a real problem it is 
more efficient and systematic to proceed as follows. Let 
us denote by X {y) the rotation by 7r/2 of all spins about 
the a;(y)-axis. As 

^-iSHyit+in+l/2)S)/2 ^ ^^tg-i<5^fy(« + ("+l/2)<5)/2;j.;^.t 

= Xe^'^"'- (*+("+l/2)5)/2 ;j.t ^ (2.7) 

it is clear that the action of e-^*^«(*+("+i/2)<5)/2 ^an be 
computed by applying to each qubit, the inverse of X 
followed by an interaction-controlled phase-shift and X. 
The prime in (2.7) indicates that Ji.j,z, ^i,2,0j ^j.z.i and 
fi^z in Hz{t + (n + 1/2)6) have to be replaced by Jij,y, 
hi,y,o, and fi^y respectively. A similar procedure 

is used to compute the action of e~"^^==(*+("+i/2)*): We 
only have to replace A" by 3^. 



By construction our algorithm to solve the TDSE (2.1) 
for spin model (2.2) is a QA itself. As a real QC 
operates on all qubits simultaneously the operation 

counts for e-**«-(*+("+l/2)5), g~»5ff„(t+(n+l/2)5)/2 g^^^ 



-i5H^(t+(n+l/2)S)/2 



0{P). 



are 0{{P + 2)), 0{{P + 2)), and 



§3. Quantum Computation 

Using the QA outlined above quantum spin systems 
containing up to 24 S=l/2 spins cau_jeasily be simu- 
lated on present-day supercomputers. EiP Here our aim 
is to use this QA-lo^simulate a recent realization of a 
2-qubit NMR Qdii-y-BEP and to execute QA's on this 
simulator. In our calculations we will take the model 
parameters corresponding to the NMR experiments of 
Refs.16,17 in which the two nuclear spins of the ^H and 
^■^C atoms in a-ijacbon-lS labeled chloroform represent 
the two qubits. Ej'EZP In the NMR set-up the molecules 
are placed in a strong static magnetic field in the +z di- 
rection. In the absence of interactions with other degrees 
of freedom this spin-1/2 system can be modeled by the 
hamiltonian 



H 



-Jl,2,zSiS2 — hi,z,oSi — /l2,z,o5'2 



(3.1) 



where /ii,z,o/27r « SOClMHz, /i2,2,o/27r « 125MHz, and 
Ji,2,z/27r « — 215Hz.t3 As the antiferromagnetic inter- 
action between the spins is much weaker than the cou- 



pling to the external field and (3.1) is a diagonal matrix 
with respect to the basis states chosen, the ground state 
of (11) is the state with the two spins up (|TT))- We 
denote this state by |00) = |t) «) |t) = ITT), i-e. the 
state with spin up corresponds to a qubit |T). A state 
of the 7V-qubit QC will be denoted by \x1X2 ■ ■ - Xn) — 
\xi) (K) 12:2) • ■ • \xn}- 

As usual it is expedient to write the TDSE for 
this NMR problem in frames of reference rotat- 
ing with the nuclear spin. Substituting |$(i)) = 
git('ii..,oSf-hh2,.,oS|)|^(-^)^ ^jjg yj^g evolution of ^'(i) in 

the absence of RF-fields is governed by the hamiltonian 
H = —Ji,2.zSiS2- This transformation has no effect on 
the expectation values of z-components of the spins but 
leads to oscillatory behavior of the x or y components, 
reflecting the fact that the spins arc rotating about the 
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z-axis. In the following it is implicitly assumed that the 
basis states of the spins refer to states in the correspond- 
ing rotating frame, even if we use the same notation for 
the basis states. 

NMR uses radio£i:equency electromagnetic pulses to 
rotate the spinsB3'E3 By tuning the frequency of the 
RF-field to the precession frequency of a particular spin, 
the power {— intensity times duration) of the applied 
pulse controls how much the spin will rotate. The axis 
of the rotation is determined by the direction of the ap- 
phed RF-field (see Refs. ||,|2|). 

§4. Grover's database search algorithm 

Finding a particular entry in an unsorted list of N 
elements is a basic problem of searching databases. In 
general this takes of the order of N operations on a con- 
ventional computer. It has been shown that a QC can 
find the item using only O (v^^ attempts.i'& 

Consider the extremely simple case of a database con- 
taining four items and functions fi{x) that upon query 
of the database return minus one for the particular item 
we are searching for and plus one otherwise. Assum- 
ing a uniform probability distribution for the item to 
be in one of the four locations, the average number 
of queries required by a conventional algorithm is 9/4. 
With Grover's QA the correct answer can be found in a 
single query (this result only holds for a database with 
4 items). Grover's algorithm for the four-item database 
can be implemented on a 2-qubit QC. 

The key ingredient of Grover's algorithm is an oper- 
ation called "inversion about the mean" that replaces 
each amplitude of the basis states in the superposition 
by two times the average amplitude minus the amplitude 
itself. This allows then for the amplification of the am- 
plitude of the basis state that represents the searched-for 
item. To see how this works it is useful to consider an 
example. Let us assume that the item to search for cor- 
responds to e.g. number 2 (/2(0) = /2(1) = /2(3) = 1 
and /2(2) = —1). Using the binary representation of in- 
tegers with the order of the bits reversed, the QC is in 
the state (up to an irrelevant phase factor as usual) 



i*) = 7T(in) + UT>-in) + iu))- 



(4.1) 



We return to the question of how to prepare th is st ate be- 
low. The operator D that inverts states like ([4.l[ ) about 
their mean reads 



D 




ITT) 

UT) 
ITi) 
lii) 



(4.2) 



The mean amplitude (i.e. the sum of all am plitudes di- 
vided by the number of amplitudes) of (4T) is 1/4 and 
we find that 



^1*) = ITi), 

i.e. the correct answer, and 
1 



:(ITT) + UT) + ITi) + Ui)), 



(4.3) 



(4.4) 



^'1*) = -^(lTT> + UT> 



lTi> + IU» = -|*>,(4.5) 



showing that (in the case of 2 qubits) the correct answer 
(i.e. the absolute value of the amplitude of ||t) equal to 
one) is obtained after 1, 4, 7, ... iterations. In general, 
for more than two qubits, more than one application of 
D is required to get the correct answer. In this sense the 
2-qubit case is somewhat special. 

The next task is to express the preparation and query 
steps in terms of elementary rotations. For illustrative 
purposes we stick to the example used above. We assume 
that initially the QC is in the state with both spins up 
(ITT))- We follow the convention used earlier in this pa- 
per, i.e. the one used in Ref. Qq) and therefore deviate 
from the notation used in Ref. ]l7|). Transforming to 
the linear superposition (4.1) is a two-step process. First 
the QC is put into the tmiform superposition state: 



\U)^W2Wi\^]) 
where 



-^(ITT) 



UT) + ITi> + Ui)), (4.6) 



V2 U 



V2 I * 



and 



"V2I1 



(4.7) 



(4.8) 



(4.9) 



represents the Walsh-Hadamard (WH) transform on 
qubit 2 which transforms ||) to i(|t) + ||))/-\/2, a clock- 
wise rotation of spin j about 7r/2 around the a;-axis and 
a anti clock-wise rotation of spin j about 7r/2 around 
the 2/-axis respectively. The inverse of a rotation Z is 
denoted by Z. The second step is to encode the infor- 
mation in the database, e.g. in the state of the 
QC. This can be accomplished by a transformation Fi 
that in the case of our example /2 [x) takes the form 





1 
-1 




ITT) 
UT) 
ITi) 
lii) 



(4.10) 



This transformation can be implemented by first letting 
the system evolve in time: 



m\v) 



(ITT) + liT) + lTi) + lii)) 



= ^(e-"/'lTT) + e+''^/'liT) 



,+i7r/4 



ITi) 



-in/i 



lii))- 



(4.11) 



For the NMR-QC based on hamiltonian (|]^) this means 
letting the system evolve in time (without applying 
pulses) for a time tq = —it/ Ji^2,z (recall Ji,2,z < 0). 

Next we apply a sequence of single-spin rotations to 
change the four phase factors such that we get the desired 
state. The two sequences YXY and YXY are particu- 
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lary useful for this purpose. We find 



\^ V V V V \^ 


o(.'= ITT) 








+ e-^"/1iT> 






p-iTr/4 




= ^-(ITT)^ 


-UT)-|TT)^ 



UT> 



lii))- (4.12) 

Combining (4.11) and (4.12) we can construct the se- 
quence transforming the uniform superposition in the 
state that corresponds to fi{x): 

Fo = YiXiYiY2X2Y2liTr), (4.13a) 

^^1 = YiXiYiY2X2Y2l{Tr), (4.13b) 

F2 = YiXiYiY2X2Y2l{Tr), (4.13c) 

F3 = YiXiYiY2X2Y2l{Tr). (4.13d) 

The remaining task is to express the process of inver- 



sion about the mean, i.e. the matrix D (see (4.2)), by 
a sequence of elementary operations. It is not difficult 
to see that D can be written as the product of a WH 
transform, a conditional phase shift P and another WH 
transform: 



D = W1W2PW1W2 
( ^ - 





V 



W1W2. (4.14) 



The same approach that was used to implement f2{x) 
also works for the conditional phase shift P (= — ^o) 
and yields 



P :^YiXiYiY2X2Y2l{T^)- 
The complete sequence Ui reads 

U, = WiW2PWiW2Fi. 



(4.15) 



(4.16) 



Each sequence Ui can be shortened by observing that in 
some cases a rotation is followed by its inverse. As there 
are various representations of the WH transform Wi that 
accomplish the same task, the sequence for e.g. i — 2 can 
be written as 

W1W2F2 = -XiXiY^X2X2Y2YiXiYiY2X2Y2l{TT) 

= -XiYiX2Y2I{t:). (4.17) 

The sequences for the other cases can be shortened as 
well, yielding 

Uo = XiYiX2Y2lin)XiYiX2Y2lin), (4.18a) 

Ui = XiYiX2Y2lin)XiYiX2Y2lin), (4.18b) 

U2 - XiYiX2Y2lin)XiYiX2Y2lin), (4.18c) 

U3 - XiYiX2Y2l{n)XiYiX2Y2l{n), (4.18d) 

where in Ui and C/2 we have dropped a physically ir- 
relevant sign. Note that the binary representation of i 



translates into the presence (1) or absence (0) in (4.18) 
of a bar on the rightmost Xi and X2. 

Table I. Specification of the EO's of a mathematically perfect 2- 
qubit QC. The execution time of each micro instruction is given 
by the second row {t/2it). The inverse of e.g. Xi is found by 
reversing the sign of hi^x,o- Model parameters omitted are zero 
for all EO's. 





Xi 


X2 


Yi 


Y2 




t/2it 


0.25 


0.25 


0.25 


0.25 


50 X 10* 


Jl,2,z 














-10-6 


hl,x,0 


+1 














h2,x,0 





-1 



















+1 



















-1 






Table III. Final state of the qubits after running the Grover's 
database search algorithm on an ideal QC (Qi,Q2, model pa- 
rameters given in table I), on a NMR-QC {Qi,Q2, model pa- 
rameters given in table II) and on the same NMR-QC {Qi,Q2, 
model parameters given in table II) using a different, but logi- 
cally equivalent, initialization sequence. 





Uo\U) 


Ui\U) 


U2\U) 


U3\U) 


Qi 


0.000 


1.000 


0.000 


1.000 


Q2 


0.000 


0.000 


1.000 


1.000 


Qi 


0.028 


0.966 


0.037 


0.955 


Q2 


0.163 


0.171 


0.836 


0.830 


Qi 


0.955 


0.041 


0.971 


0.027 


Q2 


0.031 


0.026 


0.971 


0.972 



§5. Results 

We now consider two different implementations of the 
Grover's QA described above. The first and most obvi- 
ous one is to make use of theoretically ideal unitary trans- 
formations to perform the EO's. One of the many pos- 
sible choices for the model parameters that correspond 
to this case are given in Table I. The second will be 
physical, i.e. we will use the simulator to carry out the 
NMR-QC experiment itself. The list of model parame- 
ters can be found in Table II. The values of the qubits 
(Qi,<52) are given by 



1 



{S! 



(5.1) 



The numerical values of the qubits in the final state 
as obtained by running Grover's QA on the simulator 
are summarized in Table HI. ^From the data in the first 
two rows it is evident that this QA performs as expected 
when the ideal EO's are used. In the ideal case, the final 
state (Qi, Q2) is the binary representation of the integer 
index of the "-1" item. 

The third and fourth row contain the data for the 
NMR-QC case. Using RF-pulses instead of ideal trans- 
formations to perform 7r/2 rotations leads to less certain 
answers: The final state is no longer a pure basis state 
but some linear superposition of the four basis states. 
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Table II. Specification of the elementary operations implementing a 2-qubit NMR QC, using the notation of Table 1. Note that a 
RF-pulse along the y{x) direction corresponds to a rotation about the x{y). Model parameters omitted are zero for all EO's. 





Xi 


X2 


Yi 


y2 


J-/ \ 


T j ZTT 


10 


An 
4U 


10 


An 
40 


50 X lU 




— 10~^ 


— 10~® 


— 10~^ 


— 10~® 


— 10~® 


hi,z,o 


1 


1 


1 


1 


1 


h2,z,0 


0.25 


0.25 


0.25 


0.25 


0.25 










0.05 


—0.05 





h2,x,l 








0.0125 


-0.0125 





fl,x 








1 


0.25 





h,x 








1 


0.25 







-0.05 


0.05 













-0.0125 


0.0125 











h,y 


1 


0.25 











h,y 


1 


0.25 












Indeed, using a time-dependent external pulse to rotate 
spins only yields an approximation to the simple rota- 
tions envisaged in theoretical work. This affects the ex- 
pectation values of the spin operators. What is beyond 
doubt though is that it is easy to read off the correct an- 
swer from the expectation values of the qubits. Clearly 
the simulator reproduces the experimental results and 
the QA seems to return the correct answer. 




1200 1400 1600 1800 



Fig. 1. Time evolution of the qubits Qi (solid line) and Q2 
(dashed line) obtained by executing Wi, W2 and sequence 
(4.18a) for the case of the NMR-QC. In all figures the time in- 
tervals for each operation have been rescaled to maJce them look 
equal. 



In an NMR experiment, application of each RF-pulse 

affects all spins in the sample. Although the response of 
a spin to the RF-field will only be large when this spin 
is at resonance, the state of the spins that are not in res- 
onance will also change. Thc;sc; successive; unitary trans- 
formations not necessarily commute with each other. If 



0.8 



0.6 



0.4 



0.2 



TV 




200 400 600 



800 1000 1200 1400 1600 1800 
Time 



Fig. 2. Time evolution of the qubits Qi (solid line) and Q2 
(dashed line) obtained by executing 14^2, Wi and sequence 
(4.18a) for the case of the NMR-QC. Interchanging the order 
in which the single-qubit operations Wi and W2 are applied 
changes the final state of the QC. 



they do not commute the state after these transforma- 
tions depends on the order in which the pulses have 
been applied. In fact the presence of these perturba- 
tions, although small, can have a devastating effect on 
the stability of the computation. This is illustrated by 
the data in the fifth and sixth row of Table III and by 
Figs. 1- 4. These results have been obtained by changing 
the order of preparing the two spin states. Instead of 
W1W2 we used W2W1 to initialize the QC, a permuta- 
tion that has no effect in the case of ideal EO's. ^From 
Table III it is clear that in the case of EO's implementing 
the 2-qubit NMR QC making this interchange leads to 
complete wrong results. Many of such examples can be 
constructed: The very fact that we cannot isolate one 
spin from the rest and perform operations on the for- 
mer only leads to phase errors that may (but sometimes 
don't) alter the outcome of the calculation completely. 
This QC architecture is intrinsically unstable to minor 
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0.4 



0.2 - 
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Fig. 3. Same as Fig.l but instead of sequence (4.18a) we used 
(4.18c). 



0.6 



0.4 



0.2 



200 400 600 800 1000 1200 1400 1600 1800 
Time 



Fig. 4. Same as Fig. 2 but instead of sequence (4.18a) we used 
(4.18c). 



modifications of the QA that are allowed from logical 
point of view, a much more severe problem than that 
of "decoherence" . We believe it might be interesting to 
investigate these instabilities experimentally. 

Can the phase errors discussed above be (partially) 
eliminated by some clever error-correction scheme? At 
present there is no indication they can: Any error- 
correction method requires adding extra spins to the sys- 
tem. The phase shift incured by the individual spins will 
contribute to the phase shifts of each of the many-body 
basis states and unless some magic cancelation takes 
place, the final result is unlikely to be more stable. On 
the other hand these unwanted phase factors are the re- 
sult of using RF pulses that only approximately imple- 
ment rotations about 90 degrees and may be reduced by 
using pulses that are more complicated than the sinu- 
soidal ones. Perhaps dissipation effects may also help 
to reduce the sensitivity to phase errors, a possibility 
that we are currently investigating. Another route to 
more stable operation might be to use a different set 
of EO's that more closely implements the ideal trans- 
formations. For instance, non-adiabatic transitions be- 



tween two levels driven by a periodic field display pecu- 
liar behaviorE3 and might be employed to manipulate 
the two-level systems. In this respect the single-Cooper- 
pair-boxE2P may hold some promise. In this solid-state 
device a non-adiabatic transition mechanism is used to 
let a Cooper pair tunnel between two states. Obviously 
there are many physical mechanisms to control the dy- 
namics of quantum spin systems. Exploring which of 
these mechanisms is useful for quantum computing may 
be a fertile area for future research. 
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